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Abstract 

The model quantum system of fermions in a one dimensional harmonic oscillator 
potential is investigated by a molecular dynamics method at constant temperature. 
Although in quantum mechanics the equipartition theorem cannot be used like in 
the Nose-Hoover-thermostat it is possible to couple an additional degree of freedom 
to the system which acts as a thermometer and drives the system towards the desired 
temperature via complex time steps. 
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1 Introduction and summary 



Statistical properties of finite systems are of great interest. The aim is to 
describe the behaviour of systems hke atomic clusters or atomic nuclei at 
finite temperatures and to investigate properties like the specific heat or phase 
transitions. These statistical properties are given by the partition function 
which for classical systems in the canonical ensemble reads 



r ^ r 1 

Z = jJ[ d'xi d'p, exp|--if 



where H denotes the Hamilton function, and for quantum systems 
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Z = tr(exp{-i^|) , (2) 

where H is the Hamilton operator (Throughout the article operators are un- 
derlined by a tilde and h = c = ks = l-)- For realistic systems like atomic 
clusters or nuclei where the Hamilton function or operator contains a (two- 
body) interaction it is hard or impossible to evaluate the partition function 
especially for the quantum description. 

Equations of motion for the investigated system are often much easier; either 
they are exactly known and can be integrated at least numerically as it is 
the case with the classical Hamilton's equation or they can be approximated 
with standard methods like TDHF or quantum molecular dynamics methods 
as it is the case on the quantum side. The idea then is to extract the desired 
thermodynamic quantities from the time evolution of the system. If the system 
is ergodic, ensemble averages can be replaced by time averages. 

During the last decade a huge progress has been made on the classical side of 
the problem (see for instance [1-4]). To put it into a few words, the basic idea 
is to exploit the equipartition theorem and use for instance the kinetic energy 
as a measure of the current temperature. The system is cooled or heated via 
pseudo-friction coefficients if the present temperature is too high or too low, 
respectively. The equations of motion are changed in the following way 



d ^ OH d ^ OH ^ ^ , ^ 

dt op^ dt axk 



|coc(E,.„(iv)-^r) . (4) 



The last line, equation (4), shows nicely how the system is driven towards 
the temperature T via the coupling to the pseudo-friction coefficient C- It 
was shown that the resulting distributions are those of the canonical ensemble 
(see e.g. the report [4]). In addition ergodicity can be improved by using more 
pseudo-friction coefficients and different couplings to the original system [3] . 

In quantum mechanics the problem is much more involved. No useful a priori 
relation between expectation values of observables and temperature like in the 
equipartition theorem can be exploited. Attempts have been made to derive 
a thermal non-linear Schrodinger equation which results in an ergodic wave 
function [5], but this method needs the relation between mean energy and 
temperature as an input, which is not known for most systems and therefore 
itself a matter of investigation. A nontrivial quantum extension of the Nose- 
method has been derived in ref. [6] which will be the subject of a forthcoming 
paper. 
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The idea of this article is to couple an additional degree of freedom to the 
original system, which serves as a thermometer. This idea was already suc- 
cessfully applied in nuclear physics in order to determine the caloric curve of 
nuclei and to investigate the nuclear liquid-gas phase transition [7] . The new 
aspect in this article is to drive the system via complex time steps towards 
the desired temperature. 

The total system is described by a time-dependent state \Q{t)) consisting of 
the original system | system{t) ) and the thermometer | thermometer {t) ) 

IQi't)) = I systemit) ) ® \ thermometer {t) ) . (5) 

According to the difference between the temperature Tth measured by the 
thermometer and the desired temperature T the total system is evolved by a 
complex time step dr. 

dr = dt{l - id(5) , d(5 oc {Tth - T)/Tth , \ Q{t) ) ^ | Q(t + dr) ) . (6) 

As can be inferred from (6) dp > results in coohng and < in heating 
of the system. 

In the present article the proposed method will be examined with the sim- 
ple but basic quantum system of fermions in a common one-dimensional 
harmonic oscillator potential. The advantage of this system is that all rela- 
tions of the canonical ensemble are either known analytically or can easily be 
computed numerically. 

In order to equilibrate the system in an ergodic sense a short-range two- 
body repulsion among the particles is used which does not destroy the implied 
picture of an ideal quantum gas [8] . 

The results of time averaging are compared to the canonical ensemble in terms 
of occupation numbers which in thermodynamic equilibrium characterize the 
ideal Fermi gas enclosed in a common harmonic oscillator potential completely. 
It is demonstrated that the system is indeed ergodic, i.e., the time averaged 
occupations numbers coincide with those obtained in the canonical ensemble. 



2 Method and setup 
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2. 1 Model state and equations of motion 



In many cases the time-dependent Schrodinger equation cannot be solved. 
The following time-dependent quantum variational principle (TDVP) [9] 



5 ['dt{Q{t) \ 1-T--H \Q{t)) = (7) 
Jti dt ~ 

allows to derive approximations to the time-dependent Schrodinger equation 
on the level of accuracy one needs or can afford. For the variation of the 
trial state {Q{t)\ in the complete Hilbert space the TDVP leads to the 
Schrodinger equation. For trial states given in terms of parameters Q{t) = 
iQuit)]!^ = 1, 2, ... } the TDVP leads to Euler-Lagrange equations which in 
their most general form can be written as 



EAAQ) = {Qit) I H\Q{t)) (8) 



where Cq = {Q{t) \ i-^ \ Q{t) ) and A^y{Q) is a skew symmetric matrix. The 
time evolution of the parameters then defines the time dependence of the 
many-body state | Q{t) ). If for instance | Q{t) ) is chosen as a Slater deter- 
minant of arbitrary single-particle states one ends up with time-dependent 
Hartree-Fock [10]. 

In this investigation the molecular dynamics approach of Fermionic Molecular 
Dynamics (FMD) [11] is taken where | Q(t) ) is a Slater determinant of single- 
particle Gaussian wave-packets | qi {t) ) 



(f |g;(^)) = exp|- '^a^l) J ^ , h = n + lam . (9) 

Each single-particle state is parametrized in terms of the time-dependent 
mean position mean momentum pi{t) and the complex width aiit). In the 
notation of (9) the vector hi is composed of r/, pi and a^. The time dependence 
of the spin degrees of freedom is not considered, instead it is assumed that all 
particles are identical fermions with the same spin component 1 1 ) .In order 
to perform complex time steps the complex parameters hi and ai are used. 

The approach of FMD has several advantages, in connection with this inves- 
tigation it is important to realize that for the one-body Hamiltonian 
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^ = EMO, h{l)^^ + -muY{l) (10) 

of the harmonic oscillator the solution of the approximate equations of motion 
(8) coincides with the exact solution of the Schrodinger equation. For fermions 
in a harmonic oscillator these equations of motion can be given explicitly [8] 



^bi = -imu'^aibi , ^ ai = -imu'^af + —. (11) 
at at m 



Replacing the complex hi in eq. (11) yields 



— n = — and -r.Pi^ -rnw ri . (12) 
dt m dt 

The parameters r; and pi follow the classical trajectories. Nevertheless, the 
parametrized trial state is the exact solution of the Schrodinger equation. It 
is also worth to note that for this example the equations of motion of the 
parameters are the same, regardless whether the many-body state is antisym- 
metrized, symmetrized or simply a product state. 



2.2 Time averaging 

The time averaging of an operator B is done by taking the trace with the 
statistical operator ^erg [8] 

^erg lim — ^ fUt \Q{t)){Q{t)\ (13) 

*2-»oo - ti) Jtl 

which yields 

(S7 := TV (^erg b) = ^lim ^^-1^ £ dt { Q{t) \ B \ Q{t) ) . (14) 

In general the statistical operator i?crg is a functional of the initial state 
|(5(ti)), the Hamilton operator and the equations of motion. If the er- 
godic assumption is fulfilled, the statistical operator should only depend on 



5 



2.3 Thermometer and coupling 



In order to determine the temperature of the system it is coupled to a ther- 
mometer. A thermometer is a device for which thermodynamic relations, like 
those of kinetic energy and temperature in the classical known. Here 

the thermometer is realized in terms of distinguishable particles in another 
harmonic oscillator potential. 




Fig. 1. Arrangement of fcrmion wave-packets in their oscillator potential 
(solid lines) and thermometer wave-packets in their oscillator (dashed 
lines) . 

Figure 1 sketches the setup. The fermion system consists of four Gaussian 
wave-packets (thick solid line) enclosed in their harmonic oscillator potential 
(thin solid hne). The two thermometer wave-packets (thick dashed lines) re- 
side in the other harmonic oscillator (thin dashed line). The complete system 
is represented by the state \Q{t)) 

\Q{t))^ \ Qf{t))^ \ Qth{t)) 

I Qfit) ) =4 ( I Qfiit) ) ® I QMt) ) ® I g/sW ) I qf4{t) )) (15) 

\Qth{t))= \qthiit))(S) \qth2{t)) 

where \Qf{t)) describes the fermions {A being the antisymmetrization oper- 
ator) and I Qth{t) ), which is a product state, the thermometer wave-packets. 
The time evolution according to (8) is driven by the Hamilton operator 

H = Hf + Htn + Yfj + Yf,, (16) 
where Hf is the Hamilton operator of the fermion subsystem, Hth the Hamil- 
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ton operator of the thermometer subsystem, 



(17) 



ujf = 8 MeV , ro 



i<j 




Vo = (10'...10^)a;/ , 



Vfj the short-range repulsion between all fermions needed to equilibrate the 
fermion subsystem and Vf^th the momentum dependent short-range repulsion 
between fermions and thermometer wave packets, which is given in terms of 
its matrix elements in the appendix. 

There are several conditions that need to meet in order to let the method 
work. The original system better has to be ergodic by itself and the couphng 

of the thermometer must lead to thermal equilibrium between the system 
and the thermometer. These conditions depend on the interaction which is 
needed to equilibrate the system. Some experience [8,7] and physical reasoning 
show that such interactions should be short ranged, repulsive and possibly 
momentum dependent. In addition the time evolution has to run over many 
typical cycles of the system so that enough collisions happen. I also found that 
the thermometer works best if Uth > ujf so that the thermometer has a smaller 
heat capacity and that the mass rrith < m/ so that the thermometer wave 
packets cross the region of the fermions. I chose cuth — T^ojf and rrith = m//7r^. 

If the system is ergodic the total energy of the thermometer subsystem deter- 
mines the temperature Tth 



3 Results 

3.1 Canonical ensemble of fermions in a harmonic oscillator 

Since all results of the time averaging are compared to the canonical ensemble, 
the canonical ensemble is shortly explained for N fermions in a common one- 
dimensional harmonic oscillator potential Hf. The statistical operator is given 




(18) 



by 
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5 10 15 20 5 10 15 20 

n T (MeV) 



Fig. 2. Canonical ensemble for four fermions: Occupation numbers p„ 
of the oscillator eigenstates for five temperatures (l.h.s.) and excitation 
energy as a function of temperature (r.h.s). 



5(T) = ^exp{-f} (19) 



N 



1=1 n=0 



c ' c 



then the statistical mean of an operator B is calculated as 



((S)) . :=Tr (i?(T)S) 



(20) 



In eq. (20) the subscript T indicates that the average in the canonical ensemble 
is taken at a constant temperature T. 

The most prominent physical quantities which will be investigated in the fol- 
lowing are the mean energy, which can be represented in closed form as 



((^/)) ^=^o(A^)+i:^Y 

fe=i 



cothi k— 
2T 



- 1 



(21) 



and the mean occupation probabilities p„ 



Pr, 



(22) 



where cj'^ denotes the creation operator of a fermion in the single-particle 
energy-eigenstate | n ) . The occupation probabilities are displayed in fig. 2 
on the l.h.s. for various temperatures, the r.h.s. shows the dependence of the 
mean excitation energy on the temperature. 
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5 10 15 20 " 2 4 6 5 10 15 20 

n t (fm/c) xlO* n 



Fig. 3. Occupation probabilities pn, l.h.s.: distribution sd, t = 0; middle: 

time evolution Pnit) for n = — top line, n = 3 — middle, n = 6 — 
bottom line; r.li.s.: time averaged distribution (symbols), the line shows 
the canonical result. 



3.2 Time evolution in the thermostat 



The deterministic time evolution of the system is followed over about 400 
cycles (= 60000fm/c) of the fermion system. Initially the four-fermion system 
is far from equilibrium which can be inferred from fig. 3 (l.h.s.) where the 
occupation probabilities of an initial state \Qf{t — 0)) are displayed. 

Due to the interaction between all wave packets and due to the heating or 
cooling according to the present temperature the occupation numbers are re- 
shuffled with time. This is seen in fig. 3 (middle) where the time evolution of 
the momentary occupation numbers Pn{t) is presented for n = 0, n = 3 and 
n = 6. The picture reminds of deterministic chaos. The result of time averaging 
is shown on the r.h.s. of fig. 3 (symbols) where the dependence obtained in the 
canonical ensemble is given by the sohd line (see also fig. 2). How the system 
evolves from non-equilibrium towards equilibrium can be visualized nicely by 
the occupation number of state 3 which in all three figures is given in bold. 

_i , I , I i , I : 

2 4 6 

t (fm/c) xio^ 

Fig. 4. Time evolution of actual temperatures (dashed lines) and time 
averaged temperatures (solid lines) for the thermometer system (thick 
lines) and for the fermion system (thin lines). 

Since the relation between temperature and mean energy is known for the 
fermion system (sec cq. 21), it can be checked how well the fermion tem- 
perature is reflected by the thermometer. Figure 4 shows four temperatures, 
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two for the thermometer (thick hnes) and two for the fermion system (thin 
hnes). The sohd hnes are related to the time averaged energies of the fermion 
and thermometer subsystems, respectively, where the averaging is performed 
from the beginning until the current time. The dashed lines show the momen- 
tary temperatures. One sees that the fermion system is quite excited initially 
whereas the thermometer starts being cold. One also realizes that the mo- 
mentary temperatures fluctuate strongly, which is not astonishing since both 
subsystems are very small. Taking l/\fN as a rough estimate for relative fluc- 
tuations it becomes clear that for instance in the fermion subsystem [N — 4) 
fluctuations of 50% of the mean value are quite normal. 
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Fig. 5. Time averaged occupation numbers (symbols) compared with 
those of the canonical ensemble (solid lines). 

The simulation explained in detail for T = 10 MeV is repeated for other tem- 
peratures. The resulting occupation probabilities are displayed in fig. 5. It is 
impressive that the time averaged occupation numbers match those of the 
canonical ensemble for a wide range of temperatures almost exactly. 



( ) ^ (( £n £n )) ^ Vn , (23) 

meaning that the system is ergodic, i.e., that time averages approach ensemble 
averages of the canonical ensemble. 

Nevertheless, the proposed method has also limitations. The necessary inter- 
action of all constituents prevents accurate investigations close to the ground 
state, because the system is alway excited due to the interaction. For a quan- 
tum system close to the ground state this leads already to a large temperature 
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deviation (see fig. 3). Methods to circumvent this problems are under investi- 
gation. 



A Interaction of fermions and thermometer wave packets 



The interaction between fermions (index k) and thermometer wave packets 
(index /) is modeled via the following expression 



0.01 

akR am J ^J^^S^^ 



Rl = rl + 0.2 + ^1 , ro = , V, = (10"^ fm^ . . 1 fm^) . 

Acknowledgments 



I would like to thank Prof. Klaus Barwinkel and Christian Schroder for reading 
the manuscript. 



References 



[I] S. Nose, J. Chem. Phys. 81 (1984) 511 

[2] W.G. Hoover, Phys. Rev. A31 (1985) 1685 

[3] D. Kusnezov, A. Bulgac, W. Bauer, Ann. of Phys. 204 (1990) 155 

[4] S. Nose, Prog, of Theor. Phys. Suppl. 103 (1991) 1 

[5] D. Kusnezov, Phys. Lett. A184 (1993) 50 

[6] M. Grim, E. Tosatti, Phys. Rev. Lett. 62 (1989) 2889 

[7] J. Schnack, H. Feldmeier, Phys. Lett. B409 (1997) 6 

[8] J. Schnack, H. Feldmeier, Nucl. Phys. A601 (1996) 181 

[9] P. Kramer, M. Saraceno, Lecture Notes in Physics 140, Springer, Berhn (1981) 

[10] A.K. Kerman, S.E. Koonin, Ann.Phys. 100 (1976) 332 

[II] H. Feldmeier, Nucl. Phys. A515 (1990) 147; 

H. Feldmeier, J. Schnack, Prog. Part. Nucl. Phys. 39 (1997) 393 



11 



